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I. INTRODUCTION 

Cold atoms in periodic potentials formed by stand- 
ing wave laser beams offer a test bench for a multitude 
of physics phenomena ranging from single particle band 
structure and Bloch oscillations over artificial gauge po- 
tentials to many-body transport properties and phase 
transition dynamics [IH3] . The system offers control over 
particle density and tunneling and interaction strengths, 
and read-out is accommodated by fluorescence detection 
of the atoms, either in the far field interference after re- 
lease from the lattice potential [4| or within the lattice 

urn- 

In a periodic potential, the quantum state of two atoms 
is separable in total and relative coordinates, and one has 
analytical access to states bound by attraction between 
the atoms and also to states held together by the combi- 
nation of a repulsive interaction and the band structure 
due to the lattice potential [TUTU]. A number of pub- 
lications have dealt with the separation of the center- 
of-mass and the relative motion in degenerate quantum 
gases [T0HT4] . Recently, we have [15J investigated the 
lattice system with periodic boundary conditions in the 
tight binding approximation and found that an accurate 
diagonalization of the many-body Bose-Hubbard Hamil- 
tonian leads to eigenstates which can be recognized as 
superpositions of translated replicas of a single bound 
composite many-body state. The phase factors chosen 
for this superposition govern the center-of-mass momen- 
tum of the atomic ensemble, while the relative motion of 
the atoms is accounted for by the bound composite quan- 
tum state. In [15J we verified that for sufficiently strong 
attraction, the motion within the composite object oc- 
curs on a more rapid time scale than the center-of-mass 
motion, justifying the separation of the two degrees of 
freedom for both the ground state and the lowest excited 
states of the system. 

It is the purpose of this manuscript to investigate the 
validity of a separation of the relative and center-of- 
mass motion for the problem of two attractively inter- 
acting bosons in an optical lattice similar to the Born- 
Oppenheimer approximation used in molecular chem- 
istry. In addition to the lattice we apply a confining 



harmonic potential, so that the separation of coordinates 
is not guaranteed by symmetries and conservation laws 
but has to be justified by a physical argument valid only 
in appropriate limits. 

In Sec. II, we present the Hamiltonian describing our 
system, and we derive an expression for the Hamilto- 
nian in continuous quasi-momentum space rather than 
in the discrete lattice position space. In Sec. Ill, we in- 
troduce our separation of the problem in center-of-mass 
and relative quasi-momentum coordinates, and we iden- 
tify the symmetries and boundary conditions of the states 
on the suitable reciprocal lattice. In Sec. IV, we motivate 
the Born-Oppenheimer separation of the problem, which 
leaves us with two one dimensional eigenvalue equations. 
In Sec V we present numerical solutions to the prob- 
lem, that we compare with solution of the full two-body 
Schrodinger equation. Structures in the solutions and 
spectra can be interpreted via the Born-Oppenheimer 
separation, which also offers analytical approximations 
in the different parameter limits. 



II. LATTICE HAMILTONIAN 

A. One-body Hamiltonian and Wannier states 

We consider a particle moving in a sinusoidal potential, 
so that the Hamiltonian can be written 

Aat = P 2 + ^osin 2 (^X) (1) 

where the scaled position and momentum operators have 
the dimensionless commutator 

[X,P]=L 

Since the potential is periodic with unit period, Bloch's 
theorem ensures that we can choose energy eigenstates 

with quasi-momenta q g] — 7r, it] and band indices n = 
0, 1, . . .. Another basis — the Wannier states — can be ob- 
tained as the Fourier transform over a single Brillouin 
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zone of the eigenstates 
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For sufficiently deep lattices, the Wannier states with dif- 
ferent n are localized around different lattice potential 
minima and are identical up to translation. We note 
that, within each energy band, the overlap between the 
quasi-momentum eigenstates and the Wannier states 



0\ 



Jm,n c ~ikq 



(2) 



is similar to the usual overlap between position eigen- 
states and momentum eigenstates. 

The Hamiltonian is block-diagonal in the basis of Wan- 
nier states, and the coupling of Wannier states at differ- 
ent locations is given by 



(w) m) \m at \wi n) } = -5 m , n J l 



where 



This shows that is the Fourier transform of the en- 
ergy bands as a function of q and the dispersion relations 
can be written as 



E (n: 



(n) -ikq _ 



-J, 
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2j24 n) cos(kq). 
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k=l 



For deep potentials the energy bands are relatively fiat, 
and the higher order cosine terms are suppressed. This 
justifies the tight binding approximation in which one 
retains only the nearest lattice site coupling, and in the 
following we will suppress the band index (n), and fo- 
cus on the lowest band described by the tight binding 
Hamiltonian 



HtB = ~Ji i\ W k-l) ( W k\ + \ w k+l) (W k \} 

k— — oo 

= -2Jicos(P). 

B. Harmonic confinement 



(3) 



Adding a harmonic confinement to the lattice potential 
is adequately described by adding the term kX 2 with the 
spring constant k to the Hamiltonian. For deep lattice 
potentials, the Wannier state \ wj) is well localized at X = 
j, so we make the approximation to replace X by the 
discrete quasi-position operator of the lowest band 



W 



j\Wj){Wj\ 



Introducing a rescaling of the Hamiltonian by 4Ji and 
defining k = k/AJi we end up with the Hamiltonian 



H = 



TB 



4J X 4Ji 



kW 2 . 2 cos(P) 



Similar to the usual relationship between continuous 
position and momentum operators, the discrete position 
operator W acts as a differentiation in the continuous 
quasi-momentum representation 

{%\W\a)=\^{%\a). 

which is easily derived by inserting a resolution of the 
identity in Wannier states and using the overlap for- 
mula ([2]). Therefore, we arrive at the quasi-momentum 
expression of the single particle Hamiltonian 



(%\H\a) 



" dq 2 



cos(q) 



(ip q \a) 



(4) 



At this point we make the curious observation [16j \T7\ 
that, after having restricted the Hilbert space to the low- 
est energy band and having added a quasi-harmonic con- 
finement, the Hamiltonian in momentum space Q has 
the same form as the original optical lattice Hamilto- 
nian in position space. In both cases, the Schrodinger 
equation takes the form of the Mathieu equation, but 
contrary to the case where we look for eigenstates 
with any quasi-momentum, here we will only look for pe- 
riodic eigenstates for Q, i.e. with zero "quasi-position". 



C. Interacting bosons 

In an ultra-cold gas of bosons, the interaction between 
the particles is adequately described by the two-particle 
contact interaction operator Z7i n t with the matrix ele- 
ments 

(Xi; X 2 |/7int|^3; X4) 

= gS(X 1 - X 3 )S(X 2 - X 4 )S(X 3 - X 4 ) 

for some interaction strength g. In the tight binding ap- 
proximation, the Wannier states are localized at different 
lattice sites, and one may neglect matrix elements of the 
interaction potential with Wannier product states located 
on different sites. We thus end up with the following ef- 
fective interaction operator acting on two-particle states 



Uint = G Y1 \ W ^ W j) ( W ji W 3 I 
3 



(5) 



where the strength parameter is given by 

G = g JdX \w (X)\\ 

A more rigorous treatment of the parameters of the Bose- 
Hubbard model can be found in e.g. |18j . Using the 
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relation we can calculate the matrix elements of the 
effective interaction operator in quasi-momentum space 



eff 



G 



<74 - 4i ~ 42). 



which shows that the interaction conserves the total 
quasi-momentum and is independent of its value. 

A system of two identical bosons in an optical lattice 
with harmonic confinement, which interact by the con- 
tact interaction is described by the Hamiltonian 



H 

with U = 



COs(Pi) COs(P2) 



+ U (6) 



III. RELATIVE- AND CENTER-OF-MASS 
QUASI-MOMENTA 

In the quasi-momentum representation, the cosine 
terms of Q can be written as 



cos (Pi) + cos(P 2 ) 



1^71; ^72) 



cos(gi) + cos(g 2 ) 




l^i;^ 2 ) 



I^i5^ 2 ) • 



where we have defined new operators by their action on 
quasi-momentum eigenstates, 



,iQ±/2 



The introduction of these operators suggest to re- 
parameterize the quasi-momentum basis states \^ qi ;^q 2 ) 
in terms of their sum and difference: 

4± = 42 ±gi. 

The quasi-momentum eigenstates states are defined for 
pairs of q\ and q 2 in the set 



S12 



7r; 7rJ x 



7r;7rj, 



corresponding to a diamond shaped area in the coordi- 
nate plane of q± as shown in figure [T] If we choose the 
values of in the set 

S± =] -ttjtt] x ] -27r;27r], 

then each point from £12 is represented exactly once as 
is evident from figure [T] This means that we can re- 
parametrize the quasi-momentum eigenstates as 

l4+,4-) = ^ |^(g + -g_)/2;^(g + +g_)/2) 
\^qi'^q 2 ) = y/2\Ql +42,42 ~ 4l) 
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Figure 1: (color online) Quasi-momentum of the two particles 
vs. relative and center-of-mass quasi-momentum. Left: The 
first Brillouin zone S12 in the (qi, ^2)-plane is emphasized and 
repeated in each direction. The color coding indicates the 
values of a function that is periodic in both variables with 
period 2tt and illustrates the required periodicity. The set S± 
which contains exactly one representative of each point from 
Si 2 is shown by the gray rectangle. Right: The same function 
is shown but in the (g+, g-)-coordinate system. The set S± 
is emphasized and repeated, but with a different tiling than 
for Si 2 in the left panel. 



where the front factor is chosen to preserve orthonormal- 
ity, such that we have the resolution of identity 



dq+ / d<?_ \q+,q-) (q+,q~ 

-7T J-27T 



(7) 



The corresponding discrete relative and center-of-mass 
position operators 



W± = 

act in the following way 



W 2 ± Wt 



d 



(q+,q-\W±\a) = i- — (q+,q-\a) . 
oq± 

and the interaction operator U has the following rep- 
resentation in terms of the relative and center-of-mass 
quasi-momentum states 

(q+,q-\U\a) = 7/ dq'_ (q+,q'_\a) 

J-27T 

with 7 = G/(16ttJi). 
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The two-atom Hamiltonian can now be written, 

H = 2k(WI + Wl) - cos (^^j cos I + Um ^ 

The Schrodinger equation with the Hamiltonian (|8| can 
be solved accu rately for a wide range of parameters (see 
Appendix A 1 ) . The resulting eigenenergies and the wave 
functions (A2) will be used as reference for our analysis 



by the Born-Oppenheimer separation of the motional de- 
grees of freedom which will be derived in the following 
section. 



IV. BORN-OPPENHEIMER SEPARATION 



A. Derivation 



In order to separate the relative and the center-of-mass 
motion of the system, we write the Hamiltonian in (|8| as 



H = H- 



(9) 



where H- contains all operators dealing with the relative 
motion: 



H.=2kW*-cq8 lQ± jcos 



U. 



We note that e 1( ^+/ 2 commutes with H_ and we define 
their joint eigenstates \q+,ri): 



H- \q+,n) = e n (q + ) \q+,ri) 



JQ+/2 



q+,n) 



with the following orthogonality relations 

(q+, nW+, n') = S n , n >5(q+ - q' + ). 
The states |<?+,n) can be expanded 

!?+>»> = l +2 \q'_A^\q'_)\q + , q '_). 

J-27T 



(10) 

(11) 

(12) 
(13) 



and the orthogonality relation ( 12 ) implies 

/ +2 V [A^{q_)}*A { ^\q.) = 5 n , n , 

J-27T 

Any eigenstate of the full Hamiltonian ([9| can be ex- 
panded as 



[ + \ q ' + ^cW(?;)k». 

J — IT _ 



(14) 



where the expansion coefficients C( n \q+) are found by 
applying the Hamiltonian ([9| to the expanded wave func- 
tion (14) and using (flQj) 



H |V) = [ + *dq' + C {n \q'+) Un{q' + ) + 2nW 2 + } \q' + , 



n) . 



In the g_ ^representation for the state vector, the 
eigenvalue equation takes the form of coupled differential 
equations 



EY J A^\q-)C^\q + ) 

n 

= {^(q + ) - 2K^A^\q_)C^\ q+ ). 

(15) 

The goal of the following analysis is to find an ap- 
proximation for the eigenstates, which is easier to apply 
numerically and which offers insights into their internal 
structure and dynamics. To this end, we assume that the 
states |<2 + ,n), described by q- wave functions An + \q~) 
depend only weakly on the argument q+. Eliminating 

thus the partial derivatives of An + \q~) with respect to 



q+ in the evaluation of the right hand side of (15), and 



using the orthogonality of the An + \q~) functions, we 
arrive at the following approximate equation for the ex- 
pansion coefficients 

e n (q + )C^(q + ) - 2 K ° ° 2 {q+> = ECW{q+). (16) 
oq_^_ 

This has the form of a Schrodinger equation for a single 
particle in the potential e n (</+). For each energy potential 
we can find discrete eigenenergies Effi and associated 



eigenfunctions that solve (|16|) and yield approximate 



eigenstates \ipm^) for the full Hamiltonian (|9j) 

( q+ , q -\^) = C&>( g+ )4? +) (<7-) (17) 

Note the formal similarity of this reduction of the prob- 
lem with the use of the Born-Oppenheimer approxima- 
tion in chemistry. In the latter, the wave function is 
expanded as a product of wave functions in nuclear and 
electronic coordinates, and due to the large difference in 
mass and hence in energy and time scales, the electronic 
wave functions are supposed to follow changes in the slow 
nuclear coordinates adiabatically. 

In our case, the two particles have identical masses, 
and in the absence of mutual interaction, the relative 
and center-of-mass motion occur on similar time scales, 
and the Born-Oppenheimer approximation should not be 
valid. But, as we increase the attractive interaction be- 
tween the atoms, bound states are formed, and the rel- 
ative position develops a new, faster time scale given by 
the binding energy. Our separation is carried out and 
motivated in the quasi-momentum picture, where a fur- 
ther observation may be in order: a strongly bound state 
in the relative position coordinate corresponds to a very 
extended wave function in the relative momentum, while 
the center-of-mass momentum may be well defined. This 
supports the assumption that the dominant contribu- 
tion to the second derivative in (15) stems from the q+ 

wave function Cm\q+), and hence that the derivative of 
An + \q_) with respect to q+ may be neglected. 
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Since our approximate separation of the variables 
is mathematically equivalent to the usual Born- 
Oppenheimer approximation, albeit carried out in quasi- 
momentum representation rather than position represen- 
tation, we will refer to is as "the Born-Oppenheimer ap- 
proximation" in the following. 



B. Application 

Before we apply the Born-Oppenheimer approxima- 
tion, let us consider how we expand states onto the 
center-of-mass and relative quasi- momentum eigenstates. 
Every state \qb) can be expanded in both the two-particle 
quasi-momentum basis, and in the basis of relative and 
center-of-mass quasi-momenta. 



/ dgi f_*dq 2 a(q u qx) \ip qi ; ip q2 ) 



+ 7T 



While | #1,(72) and | #+,#_) are defined for (#1,32) £ £12 
and G S±, respectively, we can look for func- 

tions defined on the entire R 2 and restrict the solution 
afterwards. In this approach, the function a is periodic 
in both variables with period 27r, and this enforces /3 to 
obey the symmetry 



/3(g + + 27r,g_±27r)=/3(<z + ,g_) 



(18) 



c.f. the tiling of R 2 with replicas of S± in the right 
panel of figure [T] Thus, a necessary — but not sufficient — 
condition is that f3 is periodic in both q + and g_ with 
periodicity 4tt. We are considering bosons and the state 
must be symmetric under the exchange of the two par- 
ticles, \-> (g + ,— g_), which implies the further 
constraint 



P(q+,q-) = fi{q+,-q-)- 



(19) 



Using these arguments on (17) we conclude that we 

are looking for solutions such that An + \q~) is even and 
periodic in g_ with period 47r, and such that the prod- 
uct of Cm\q+) and An + \q_) is periodic in q+ with the 
same period. Furthermore, the product must satisfy the 
relation (18). 



1. The first Born-Oppenheimer equation 

To apply the Born-Oppenheimer approximation, we 
must first find the eigenstates of H_ and their eigenen- 
ergies, and using the formal expansion of the states (13), 
the eigenvalue equation (10) leads to the equation 



2k* -F( q+ ) <*»(«=■) 



dq 



where F(q+) = cos (%r). For each value of q+ : this equa- 
tion has the form of a Schrodinger equation with argu- 
ment <?_, and with a periodic cos(^) potential with am- 
plitude F(q+) and a non-local potential with strength 7. 
Solutions which are periodic in q_ with period 47r are 
readily found by Fourier expansion of A^ + \q_) (see Ap- 



+2tt 



+ 7/ &q'_A^\q'_) (20) 



pendix A 2 ), and these solutions can be chosen to be real- 



valued just like the zero quasi-momentum eigenstates for 
cosine potentials in position space. 

The front factor F(q+) of the cosine potential is itself 
a cosine function of q+ leading to two observations: 



1. F(q+) is an even function of q+ so Eq. (20) is un- 
altered under the transformation q+ ^ —q+. Thus 
the solutions must be identical up to a complex fac- 
tor, and since they are real-valued we can choose 
the solutions as 



(21) 



We could not have chosen a minus sign, since this 
would have made An + ^ vanish for q+ = 0. 

2. F(q+) changes to values of opposite sign when q+ 
is increased by an amount of 2ir and the cosine 



potential cos(g_/2) in (20) is effectively translated 
by half a period. For this translated potential the 
eigenvalues are the same, while the eigenfunctions 
are translated and scaled 

tn(q+) =e n (g + + 27r) (22) 
A^\q-) = f n ^+ +2 *>(<z_ ± 2tt). (23) 

The factor £ n may take the values ±1 since 
An + \q-) is real- valued for all values of q+ and <?_. 



Applying the relations (21) and (23) for q+ = —ir we get 
the relation 



4 +OT) (9-)=£n4 +OT) (9-±27r) 



(24) 



so we can determine £ n from the translational symmetries 



2. The second Born-Oppenheimer equation 



Solving Eq. (20) yields the potential e n (g + ) which is 
periodic with period 27r, and we are looking for func- 
tions Cm\q+) that are periodic in q+ with period 47r. 
Therefore, Bloch's theorem tells us that we can choose a 
complete set of solutions as 



Cl B) («4 



^/ 2 D%\q + ) 



(25) 



) where D$ is periodic with periodic 27r, and S n = 0,1. 
For S n = the solution Cm\q+) is thus periodic with pe- 
riod 27r, whereas for S n — 1, it is antiperiodic. We require 
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that the product of Cm\q+) and An +) (q-) satisfies the 
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Figure 2: (color online) Energies and eigenfunctions found by solving the two Born-Oppenheimer equations for k = 0.5 and 
7 = —0.5. Left panel: The six lowest potential curves e n (q+) found from the first Born-Oppenheimer equation. Upper panels: 

Magnification of four of the potential curves in the left panel. Lower panels: Eigenfunctions A^ + \q-) for the first Born- 
Oppenheimer equation shown for all values of q+ for the corresponding n-values. In the upper panels are shown (horizontal 
dashed blue/red lines) the two lowest energies for m — 0, 1 found from solving the second Born-Oppenheimer equation in 
the potential e n (g+) and the corresponding wave functions (solid blue/red lines). 



symmetry (18), and if we combine this with (23), we get 
the relation 



C£\q + )A%+\q_) = t n C£\q+ + 2n)A^\q_) 



from which we conclude that Cm\q+) must fulfill the 
symmetry 



CW( 9+ + 2^) = ^W(? + ). 



Comparing to (25) we see that for £ n = —1 we must 
choose 5 n = 1 and for £ n = +1, we must use S = 0. We 
can solve {w\ by Fourier expansions of and e n (see 
Appendix |A 3|) . 



V. BORN-OPPENHEIMER SOLUTIONS 
A. Wave functions 

When solving the first Born-Oppenheimer equa- 
tion ( |20| ) we find eigenvalues e n (g + ) and eigenfunctions 
(q_) for each value of q+ g] — 7r, +7r]. In the leftmost 
panel in figure [2] is shown the six lowest potential curves 
e n (q+). The lowest potential curve is well separated from 
the higher ones which lie closer. Each of the potential 
curves has an energy variation which is typically small 
compared to the energy distance between the bands, and 
in the upper panels, a magnified view of the curves are 
shown. In the lower panels, eigenfunctions An + \q_) of 
the first Born-Oppenheimer equation are shown for four 
different values of n. 

The second Born-Oppenheimer equation uses the en- 
ergies €n(q+) as potential functions in a Schrodinger 
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Figure 3: (color online) Upper panels: Quasi-momentum wave functions found by the Born-Oppenheimer approximation for 
k = 0.5 and 7 = —0.5. Columns 1-2 show the two lowest eigenstates in the lowest potential curve eo (<?+), corresponding to the 
n — column in figure [2] Columns 3-4 correspond to the n — 1 column in figure [2] and columns 5-6 correspond to the n — 2 
column in figure [2] Lower panels: The corresponding exact solutions. 



like equation, and each of the upper panels in Figure [2] 
shows the energy levels of the two lowest eigenstates 
(m = 0, 1) in these potentials along with their eigen- 
functions Cm\q+)- As we saw in the previous section, 
the function Cm\q+) should be chosen periodic or anti- 
periodic depending on the value of £ n . By studying the 
behavior of An + \q~) at q+ = ±7r one can see if £ n 
is +1 or —1 depending on whether the wave function 
A^ z7T \q_) changes sign when translated by tt or not. 

For n = 0, 2 the solutions Cm\q+) to the second Born- 
Oppenheimer equation must be periodic with period 27r, 
while for n = 1, 5 the solutions Cm\q+) must be chosen 
cm^periodic. 



Total Born-Oppenheimer solutions to the two-atom 
Hamiltonian are shown in Fig. |3j where panels (Al-2) 
correspond to the approximate solutions from the n = 
case of figure [2J panels ( A3-4) correspond to the n = 1 
case, and panels (A5-6) correspond to the n = 2 case. In 
the lower panels of figure |3j the corresponding exact two- 
atom eigenstates are shown. There is a good agreement 
between the exact and approximate solutions, especially 
for the low excitations of the lowest bands. 



B. Energies 

In figure [4] both the exact and the approximative en- 
ergies are shown for fixed n as functions of the scaled 
interaction strength 7. Except in the region where 7 is 
numerically small, there is reasonable agreement between 
the exact and the approximated energy levels. For nega- 
tive 7 there is a clear grouping of the energy levels in two 
groups: Those that are nearly constant as a function of 7 
and those that depend linearly on 7. Comparing to the 
approximate energies found by the Born-Oppenheimer 
approximation we see that the linear dependence comes 
from the fact that the position of the lowest potential 
curve varies linearly with 7, as we will see in the follow- 
ing. 



C. Approximate solution of the first 
Born-Oppenheimer equation 

To understand the behavior of the energy spectra, we 
start by analyzing the system in the limit where at least 
one of the two coefficients 7, n is (numerically) much 
larger than unity, so that we can find analytical approx- 
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Figure 4: (color online) Exact and approximate energies as 
a function of 7 for k = 1. Exact energies are plotted by the 
black dashed lines, and for n = 0,1,2,3 the energies E^ 
found from the Born-Oppenheimer equation are plotted in 
four different colors. 



imations. This limit enables us to treat the term 

-F( 9+ )cos(^) (26) 



in the first Born-Oppenheimer equation (20) as a per- 



turbation. When we neglect this term, we can choose a 
complete set of eigenfunctions as plane waves with wave 
number k/2 for k = 0, ±1, ±2, . . . and with energies 



(27) 



but only the even linear combinations are physically rel- 
evant. Note that the k does in general not coincide 
with the excitation number n as used in the first Born- 
Oppenheimer equation, where the energy curves e n were 
sorted by energy. 

The term Air^y contributes only for k = since all 
other plane waves integrate to zero in the second line 
of Eq. (20). Even when the omission of (26) is not valid, 



the integral still becomes substantial if An + \q~) has 
no nodes, whereas it is suppressed when there are sign 
changes in An + \q~)- 



In figure[4]we notice some discontinuities in the approx- 
imate energies, which can be explained in the following 
way. Due to the linear dependence of the energy for the 
k = plane wave, its energy becomes degenerate with 
the higher levels, when 7 varies. More precisely, eo will 
cross 6/c at the 7- value 



Ik 



8tt " 



(28) 



Without the symmetry requirement ( |18| ) we could find 
two families of solutions to the second Born-Oppenheimer 
equation for each potential curve e n (q+). Depending on 

the symmetries of the solution An + \q~) discussed in 
Sec. |IV B| we can only choose one of these families, and 
at each side of the energy crossing (28), we must discard 



one or the other and thus obtain a discontinuous energy 
dependence. 



Now, we turn to the term (26), the effect of which we 



will approximate using non-degenerate perturbation the- 
ory. Due to the orthogonality between the cosine func- 
tions, there are no first-order corrections. The second 
order corrections, on the other hand, give contributions 
of the form 

Ae fc (<2+) = a k F(q+) 2 

where the amplitude a k can be calculated (see Ap- 
pendix Ib]) 



K — 8-7T7 
1 

2k— 16-7T7 
1 

k «(4/c 2 -l) 



k = 0, 
k = ±1, 
otherwise. 



(29) 



This gives the perturbative approximation to the poten- 
tial curves 



~e k (q + ) + A£ k (q+) = ^k 2 + 4tt7<J m + y (1 



cos(<2 + )) 



Due to the term 47r7 in the expression for eo (<?+), the 
energies of the eigenfunctions in this potential change 
linearly with 7. For k > 1 the position of e k {q + ) depends 
less strongly on 7 and the eigenstates in these potentials 
have almost constant energy. 



D. Approximate solution of the second 
Born-Oppenheimer equation 

To analyze in more detail how the eigenenergies 
are distributed we must take a closer look at the sec- 
ond Born-Oppenheimer equation which has the form of 
a Schrodinger equation for a particle of mass h 2 /An in 
the potential e n (g+). When the above perturbative treat- 
ment is valid, this potential is a cosine with amplitude 
\ak \ /2, so in order to estimate the eigenstates and ener- 
gies, we must compare n and \a k \. In the limit where we 
can neglect the ^-dependence of the solutions 
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I 



I 




Figure 5: (color online) Exact and approximate energies as 
a function of k for 7 = —10. The black dashed curves show 
the exact energies E n , while the solid red curves show the 
approximate energies E$ found from the lowest potential- 
curve in the second Born-Oppenheimer equation. The green 
curve shows the position of the maximum of the lowest poten- 
tial curve eo(<?+) within the perturbative approximation. The 
exact ground state energy Eo, which varies with k, has been 
subtracted from all energies, and afterward, the energies is 
scaled by the energy difference between the two lowest exact 
energy levels E\ — Eq. 



can be well approximated by plane waves e imq+ / V2n 
with "box potential' -energies 



fc,0 



2 turn 



which depend quadratically on m. In the opposite limit 
where €k(q+) is a deep potential in (16), we can ap- 



proximate the cosine potential by a quadratic expan- 
sion around its minimum. The resulting equation is a 
Schrodinger equation for a particle in a harmonic oscil- 
lator of frequency 



dk\ 



For the lower part of the energy spectrum, the solutions 
are then well approximated by the usual harmonic os- 
cillator eigenstate wave functions and the energies are 
equidistantly spaced with spacing HuJk- 



• 47T7<Sfc > o • 



Ok 

2 



Q>k\ 



Figure [5] illustrates the transition between the "particle 
in a box" and the "harmonic oscillator" regimes by show- 
ing the exact and approximate energies E$ as functions 
of n for fixed negative 7. Since the harmonic oscilla- 



deep, it requires that k <C |ao| = |ft — 8^7 1~ , so to cap- 
ture the whole transition, the ft-axis is logarithmic. The 
energies are plotted after subtracting the ground state 
energy Eo and scaling by the energy difference E\ — Eq 
between the first excited state and the ground state. For 
k <C 1 the harmonic oscillator spectrum is then revealed 
as levels with unit spacing. For k — » 1, on the other 
hand, the curves become constant at 1,4,9, .. . showing 
the quadratic dependence on m (see also [19]). We note 
that there is a perfect agreement between the exact and 
approximate energies shown in the figure. In the tran- 
sition from the harmonic oscillator regime to the "par- 
ticle in a box" regime, the energy levels group in pairs, 
which have the following explanation: For a deep po- 
tential curve e n there is a significant energy difference 
between the first excited even and odd states, but when 
the potential curve is nearly constant, then even and odd 
solutions with a given wave number has almost the same 
energy. 

No matter how deep the potential curve €k(q+) is, the 
harmonic approximation is not perfect, and above some 
energy the spectrum is ill-described by a harmonic os- 
cillator spectrum. A simple estimate suggests that the 
description is good for eigenstates whose energies lie be- 
low the maximum of the potential curve, which is ap- 
proximated by the unperturbed energies (27) plus a term 
depending on the sign of a& 



tk(7,K>) = -k 2 



■ 47T7<5fc,o 



Q>k + K| 



tor approximation is valid when the potential in (16) is 



In figure [5] this (solid green) curve is shown for k = 
and agrees systematically with the border where the har- 
monic oscillator energy spectrum is significantly altered. 



VI. CONCLUSION 

In the present paper we have considered two identi- 
cal bosons on an infinite, discrete lattice with an addi- 
tional harmonic confinement. In the tight binding ap- 
proximation, the single particle physics in terms of quasi- 
momenta is described by the same equation as a sin- 
gle particle in a continuous cosine potential — namely the 
Mathieu equation. Adding a contact interaction yields 
a Hamiltonian which does not separate in relative and 
center-of-mass coordinates, even though the two-body in- 
teraction problem separates in both a homogeneous dis- 
crete lattice Hamiltonian and in a continuous harmonic 
oscillator. 

By formulating the problem in quasi-momentum repre- 
sentation we can make an approximation which is mathe- 
matically equivalent to the usual Born-Oppenheimer ap- 
proximation performed in position space in molecular 
physics: We thus find approximate solutions by first solv- 
ing an equation for the relative quasi-momentum wave 
function that depends parametrically on the center-of- 
mass quasi-momentum. This yields potential curves for 
a Schrodinger equation for the center-of-mass coordinate, 
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which is readily solved. Contrary to the usual Born- 
Oppenheimer approximation used to separate slow nu- 
clear and fast electronic motion in molecules, in our sys- 
tem we have a tunable adiabaticity parameter, namely 
the strength of the inter-particle interaction. 

In the solution of both the first and second Born- 
Oppenheimer equations we can identify the excitation 
degrees of freedom in the system. This provides physi- 
cally motivated quantum numbers valid also for the exact 
eigenstates together with rules for which quantum num- 
bers are allowed by symmetry considerations. 

Finally, from the good agreement between the exact 
and approximate solutions we conclude that the Born- 
Oppenheimer approximation is well justified when the 
energy scales for the relative and the center-of-mass mo- 
tion of the two-particle quantum state are well-separated. 
We imagine that a similar separation may be useful for 
approximate first principle calculations on many other 
cold atom systems, e.g., with more particles and possibly 
with mixtures of different species. 
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Appendix A: Solving the equations numerically 

1. Solving the non-approximated equation 

To solve the two-atom Schrodinger equation in the 
tight binding approximation we expand the state as 



\a) = ^F jk \wj-w k ) . 

3,k 



(Al) 



The stationary Schrodinger equation with the Hamilto- 
nian ([8j yields the equation for the expansion coefficients 

EF j)k = 47r-yS jik F jik + n(j 2 + k 2 )F jk 



The original Hamiltonian is invariant under parity inver- 
sion of both particles so we can find a complete set of 
solutions of even and odd wave- functions. In terms of 



the expansion ( Al ) this means that we can find solutions 
where 

F-j,-k = pFj,k 

where p can assume the values ±1. In addition, since we 
are dealing with two identical bosons, only symmetrized 
wave functions are physically meaningful, with implies 
that we have the symmetry 



For numerical purposes we enforce these requirements 
by hand in the following way. Instead of looking at all 
pairs (j, fc) E Z 2 , we restrict our attention to those in the 
subset 

T = {(j,k)ez 2 | |fc| < j< j max } cZ 2 . 

for some manually chosen j max . Using the symmetries we 
reformulate the recurrence equation such that it only in- 
volves coefficients from T. The equation can be expressed 
as a matrix eigenvalue equation which is amenable to 
standard numerical diagonalization routines. When all 
coefficients have been found — and properly normalized — 
the wave function in relative and center-of-mass quasi- 
momenta is given by 



(q+,q-\<x) = ^2 F jk (q+,q-\wj]W k ) 



(k-j)q-/2 

(A2) 



2. Solving the first Born-Oppenheimer equation 



The solutions of (20) are functions An + \q~) which are 
periodic in q_ with period 47r. Therefore, for each value 
of q+ we make the expansion 



(q+)\jq-/2 



and obtain the tridiagonal recurrence relation, 



(A3) 



a 



(fi 2 + 4^,0 -<*(?+)) 4 9 n +) - (A4) 



To accommodate the bosonic nature of the particles, we 



only look for even solutions to (20), so we only need to 
consider terms ol^^ with j > 0, and for j = we use 

F(q + )a{ q : ) = [^l-Uq + )}^t ) - ( A5 ) 



By expressing the recurrence relation as a matrix eigen- 
value equation, this can be truncated and solved with 
good accuracy. 



3. Solving the second Born-Oppenheimer equation 

We solve the second Born-Oppenheimer equation using 
the results from the first Born-Oppenheimer equation. 
First, coefficients in the expansion 



Fj,k = Fi 



k,j- 



(A6) 
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are determined by a discrete Fourier transformation. We 
then use the expansion 

c£>fao = iUE^ (,+ * ),+ - 



in (16) together with the expansion ( A6) which yields the 



following equation for the 7-coefficients 



Using the orthogonality of the cosine functions we see 
that only coefficients with neighboring values of I and k 
are coupled 



1 lk — ~ [Ok-l+1 + Gfc-J-lJ- 



£/W + 2*e(l 



7; 



(A7) The resulting perturbative corrections then take the form 



Since the potential energy curves e n (#+) are even func- 
tions, the solutions can be chosen to be either even or 
odd, and the coefficients then fulfill 7 i m ' n = ±7™' n . It 
suffices to only consider coefficients with m > and solve 
the recurrence equations. 



Appendix B: Calculation of perturbation terms 

The second order perturbation terms for the poten- 
tial curves Ck(q+) iare found by calculating the matrix 
elements of the term (26) between pairs of unperturbed 



eigenfunctions which are plane waves: 



_ F(q+) 
An 



/ C 

J-2ir 



i(fc-Qg_/2 



cos 



(y)- 



Ai k (q+) = J2 J 



\h k Y 



l^k 



Ck — £r 



a k F(q + f 



where the amplitude of the oscillation is 



[5k-i+i + ^k-i-iY 



l^k 



(k* - P) + 167T7(4,o - M ' 



Here we can distinguish between the three cases k = 0, 
k = d=l and |fc| > 2, where we get the results summarized 
in Eq. (1291). 
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